/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */



#ifndef _OS_distribution_
#define _OS_distribution_

#include <src/data/Operators.h>
#include <src/miscellaneous/OS_basic.h>
#include <src/miscellaneous/Function_Utilities.h>
#include <boost/math/distributions.hpp> // we follow the boost math library

#define OS_normal									boost::math::normal_distribution
#define OS_non_central_chi_squared					boost::math::non_central_chi_squared_distribution

template<class Real_Type = OS_double>
struct base_normal_distribution : 
	public OS_normal<Real_Type> ,
	public operator_base<Real_Type,Real_Type>
{
	typedef OS_normal<Real_Type> Base;
	base_normal_distribution(const Base& from) : Base(from)
	{};
	base_normal_distribution(argument_type average = 0, argument_type sigma = 0.05) : Base(average, sigma)
	{};
	Real_Type operator()(Real_Type x) const {
		return boost::math::pdf((Base&)(*this), x);
	};
};



template<class Real_Type = OS_double>
struct OS_normal_distribution : public template_function_utility <base_normal_distribution<Real_Type> > 
{
	typedef template_function_utility Base;
	OS_normal_distribution(const base_normal_distribution<Real_Type> & from) : Base(from)
	{};
	OS_normal_distribution(argument_type average = 0, argument_type sigma = 0.05) : Base(average, sigma)
	{};
};

/* This function is the increasing rearangement of the bounded variation density of the normal distribution.
	
	Let the normal distribution $\rho(x) = \frac{1}{\sqrt{2\pi}} e^{\frac{-x^2}{2}}.
	Let denote $\rho_{\sigma}^{\bar{\rho}}(x) = \frac{1}{\sigma}\rho( \frac{x-\bar{\rho}}{\sigma})$ the normal distribution of variance $\sigma$ and average $\bar{\rho}$.
	Let consider the positive function (bounded variation density)
	\begin{equation}
		u(x) = \frac{ |\partial_x \rho_{\sigma}^{\bar{\rho}}(x)| }{\int_{-\infty}^{+\infty} |\partial_x \rho_{\sigma}^{\bar{\rho}}(s)|ds}
	\end{equation}
	Then this distribution is the increasing rearangement of the previos distribution, that is a function $S_{\sigma}^{\bar{\rho}}(y)$ from [0,1] into $\R$ such that
	\begin{equation}
		\partial_y S_{\sigma}^{\bar{\rho}} u(S_{\sigma}^{\bar{\rho}}(y)) = 1
	\end{equation}
	equivalently, $S_{\sigma}^{\bar{\rho}}$ is the inverse function of
	\begin{equation}
		\frac{\int_{-\infty}^x |\partial_x \rho_{\sigma}^{\bar{\rho}}(s)|ds}{\int_{-\infty}^{+\infty} |\partial_x \rho_{\sigma}^{\bar{\rho}}(s)|ds}
	\end{equation}
	Computations: Let $S=S_{1}^{0}$. The dependance is
	\begin{equation}
		S_{\sigma}^{\bar{\rho}} = \sigma S + \bar{\rho}
	\end{equation}
	With elementary computations, we get
	\begin{equation}
		S(y) = \{-\sqrt{2 \ln (2y) }, y < 1/2 ; +\sqrt{2 \ln (2(y-1)) }, y >1/2   \}
	\end{equation}
*/
template<class Real_Type = OS_double>
struct base_inverse_normal_distribution : 
	public operator_base<Real_Type,Real_Type>
{
	base_inverse_normal_distribution(argument_type average = 0, argument_type sigma = 1.) : average_(average), sigma_(sigma)
	{
		OS_DYNAMIC_CHECK(sigma_>=0,"incoherent");
	};
	Real_Type operator()(Real_Type y) const {
		if (y<0.5)
			return -sigma_*std::sqrt(-2.*std::log(2.*y)) + average_;
		return sigma_*std::sqrt(-2.*std::log(2.*(1.-y))) + average_;
	};
	Real_Type average_;
	Real_Type sigma_;
};



template<class Real_Type = OS_double>
struct OS_inverse_normal_distribution : public template_function_utility <base_inverse_normal_distribution<Real_Type> > 
{
	typedef template_function_utility Base;
	OS_inverse_normal_distribution(const base_inverse_normal_distribution<Real_Type> & from) : Base(from)
	{};
	OS_inverse_normal_distribution(argument_type average = 0, argument_type sigma = 1.) : Base(average, sigma)
	{};
};


template< class Real_type = OS_double>
struct OS_base_cumulative_normal_distribution : public derivable<Real_type>
{
  public:
  typedef OS_base_cumulative_normal_distribution type;
	type(const type & from) : normal_(from.normal_) {}
	type(Real_type average = 0.0, Real_type sigma   = 1.0) : normal_(average,sigma) {}
    // function
	Real_type operator()(Real_type x) const {
		return boost::math::cdf(normal_, x);
	};
	Real_type derivative(OS_double x) const {
		return normal_(x);
	};
	virtual smart_ptr_function get_derivative() const {
		return smart_ptr_function(new OS_normal_distribution<Real_type>(normal_));
	};
	OS_normal<Real_type> normal_;
};

template<class Real_Type = OS_double>
struct OS_cumulative_normal_distribution : public template_function_utility <OS_base_cumulative_normal_distribution<Real_Type> > 
{
	typedef template_function_utility Base;
	OS_cumulative_normal_distribution(argument_type average = 0, argument_type sigma = 0.05) : Base(average, sigma)
	{};
};

template< class Real_type = OS_double>
struct OS_base_inverse_cumulative_normal_distribution : public derivable<Real_type>
{
  public:
	  OS_base_inverse_cumulative_normal_distribution(Real_type average = 0.0, Real_type sigma   = 1.0) : normal_(average,sigma) {}
    // function
	Real_type operator()(Real_type x) const {
		return boost::math::quantile(normal_, x);
	};
	Real_type derivative(OS_double x) const {
		return 1./normal_(this::oprator(x));
	};
	Real_type average() const {
		return boost::math::mean(normal_);
	};
	Real_type sigma() const {
		return boost::math::variance(normal_);
	};
  private:
	 OS_normal<Real_type> normal_;
};

template<class Real_Type = OS_double>
struct OS_inverse_cumulative_normal_distribution : public template_function_utility <OS_base_inverse_cumulative_normal_distribution<Real_Type> > 
{
	typedef template_function_utility Base;
	typedef OS_inverse_cumulative_normal_distribution type;
	type(const type & type) : Base(type)
	{};
	type(argument_type average = 0, argument_type sigma = 0.05) : Base(average, sigma)
	{};
};


template<class Real_Type = OS_double>
struct base_OS_non_central_chi_squared : 
	public OS_non_central_chi_squared<Real_Type> ,
	public operator_base<Real_Type,Real_Type>
{
	typedef OS_non_central_chi_squared<Real_Type> Base;
	base_OS_non_central_chi_squared(argument_type delta = 1., argument_type lambda = 0.0) : Base(delta, lambda)
	{};
	base_OS_non_central_chi_squared(const Base& from) : Base(from)
	{};
	Real_Type operator()(Real_Type x) const {
		return boost::math::pdf((Base&)(*this), x);
	};
};

template<class Real_Type = OS_double>
struct OS_non_central_chi_squared_distribution : public template_function_utility <base_OS_non_central_chi_squared<Real_Type> > 
{
	typedef template_function_utility Base;
	OS_non_central_chi_squared_distribution(argument_type delta = 1., argument_type lambda = 0.0) : Base(delta, lambda)
	{};
};


template< class Real_type = OS_double>
struct OS_base_cumulative_non_central_chi_squared_distribution : public derivable<Real_type>
{
  public:
	OS_base_cumulative_non_central_chi_squared_distribution(argument_type delta = 1., argument_type lambda = 0.0) : non_central_chi_squared_(delta,lambda) {}
    // function
	Real_type operator()(Real_type x) const {
		return boost::math::cdf(non_central_chi_squared_, x);
	};
	Real_type derivative(OS_double x) const {
		return non_central_chi_squared_(x);
	};

	virtual smart_ptr_function get_derivative() const {
		return smart_ptr_function(new base_OS_non_central_chi_squared<Real_type>(non_central_chi_squared_));
	};
  private:
	 OS_non_central_chi_squared<Real_type> non_central_chi_squared_;
};

template<class Real_Type = OS_double>
struct OS_cumulative_non_central_chi_squared_distribution : public template_function_utility <OS_base_cumulative_non_central_chi_squared_distribution<Real_Type> > 
{
	typedef template_function_utility Base;
	OS_cumulative_non_central_chi_squared_distribution(argument_type delta = 1., argument_type lambda = 0.0) : Base(delta, lambda)
	{};
};

template< class Real_type = OS_double>
struct OS_base_inverse_cumulative_non_central_chi_squared_distribution : public derivable<Real_type>
{
  public:
	OS_base_inverse_cumulative_non_central_chi_squared_distribution(argument_type delta = 1., argument_type lambda = 0.0) : non_central_chi_squared_(delta,lambda) {}
    // function
	Real_type operator()(Real_type x) const {
		return boost::math::quantile(non_central_chi_squared_, x);
	};
	Real_type derivative(OS_double x) const {
		return 1./non_central_chi_squared_(this::oprator(x));
	};
	Real_type average() const {
		return boost::math::mean(non_central_chi_squared_);
	};
	Real_type sigma() const {
		return boost::math::variance(non_central_chi_squared_);
	};
  private:
	 OS_non_central_chi_squared<Real_type> non_central_chi_squared_;
};

template<class Real_Type = OS_double>
struct OS_inverse_cumulative_non_central_chi_squared_distribution : public template_function_utility <OS_base_inverse_cumulative_non_central_chi_squared_distribution<Real_Type> > 
{
	typedef template_function_utility Base;
	typedef OS_inverse_cumulative_non_central_chi_squared_distribution type;
	type(const type & type) : Base(type)
	{};
	type(argument_type delta = 1., argument_type lambda = 0.0) : Base(delta, lambda)
	{};
};




#endif
